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Abstract 

We present a novel study on the problem of constructing mass models for the 
Milky Way, concentrating on features regarding the dark matter halo component. 
We have considered a variegated sample of dynamical observables for the Galaxy, 
including several results which have appeared recently, and studied a 7- or 8- 
dimensional parameter space - defining the Galaxy model - by implementing a 
Bayesian approach to the parameter estimation based on a Markov Chain Monte 
Carlo method. The main result of this analysis is a novel determination of the 
local dark matter halo density which, assuming spherical symmetry and either an 
Einasto or an NFW density profile is found to be around 0.39 GeV cm~ 3 with a l-cr 
error bar of about 7%; more precisely we find a pdm{Rq) = 0.385±0.027 GeV cm" 3 
for the Einasto profile and p D Ai{Ro) = 0.389 ±0.025 GeV cm" 3 for the NFW. This 
is in contrast to the standard assumption that pdm{Rq) is about 0.3 GeV cm~ 3 
with an uncertainty of a factor of 2 to 3. A very precise determination of the 
local halo density is very important for interpreting direct dark matter detection 
experiments. Indeed the results we produced, together with the recent accurate 
determination of the local circular velocity, should be very useful to considerably 
narrow astrophysical uncertainties on direct dark matter detection. 
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1 Introduction 



There is by now incontrovertible evidence that dark matter is the building block of all 
structures in the Universe. This statement is motivated by the successes of the currently 
favored model for structure formation, the ACDM model, in accounting for a variety of 
complementary cosmological observations, as well as by dynamical observations per- 
formed on the largest structures we see in the Universe, namely galaxy clusters, down to 
low surface brightness and dwarf galaxies. Within the picture one obtains when merging 
such wealth of information, dark matter behaves like a dissipation-less and collision-less 
fluid and it is cold, i.e. the free-streaming scale associated to dark matter is negligi- 
bly small (more precisely, there are rather tight constraints on the violation of each of 
these properties). It follows that none of the particles in the Standard Model of Particle 
Physics is suitable to play the role of dark matter. Although there are great expectation 
from the LHC to shed light on physics beyond the Standard Model and an eventual 
embedding of a dark matter candidate within it, a clean handle on the dark matter 
puzzle will come only from direct or indirect detection signal of dark matter particles 
within dark matter halos. 

There are classes of dark matter candidates for which it is indeed feasible to extract a 
direct detection signal, namely a signal of the interaction of a dark matter particle from 
the local population of the Milky Way dark matter halo within a laboratory detector: 
the most extensively studied models are those of Weakly Interacting Massive Particles 
(WIMPs) and axions. In both cases the signals scale linearly with the local number 
density of dark matter particles (and depend also on their local velocity distribution). 
It is often quoted that such density is merely known within a factor of 2 or 3, an 
uncertainty which limits significantly the possibility of cross-correlating, e.g., for a given 
super symmetric model, the production cross section at the LHC with the signal induced 
by the scattering rate of the lightest supersymmetric particle (dark matter candidate) 
on a target nucleus in an underground detector. The problem of the not very accurate 
determination of the local halo density goes back to seminal papers in the field in the 
eighties and nineties, such as the analysis by Gates, Gyuk & Turner [1,2], which are 
however based on compilations of dynamical constraints for the Galaxy that are by now 
outdated. In this respect, there has been significant improvement in the last few years, 
and it is reasonable to expect that this would have an impact on the inferred local halo 
density. 

There are very strong indications that spiral galaxies are embedded into massive 
dark halos. These include the flatness of rotation curves up to their largest measurable 
radii (see, e.g., [3]), the dynamics of the satellites (see, e.g., [4]), weak lensing results 
(see, e.g., [5]), and warping and flaring of the gas layers (see, e.g., [6,7]). It is not, 
however, straightforward to derive detailed models for dark halos from such dynamical 
observations. For external galaxies, one options is make a compilation of data for dif- 
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ferent objects, extrapolate trends to be correlated to physical properties of the object 
and constrain the mean functional form for the dark matter profile, see, e.g., the recent 
results in [8,9]. The focus here is however on a single object, the Milky Way, and on 
addressing whether, starting from a rather generic model for the Galaxy, a careful cross 
correlation of all currently available data is sufficiently powerful to provide a relevant 
constraint on quantities like the local halo density. 

Several authors have considered the problem of building mass models for the Galaxy 
(an incomplete list of references includes [2, 10-16]). The standard approach is to 
describe the Galaxy as a superposition of density profile components, and to discriminate 
among viable ansatz by comparing with available dynamical constraints. We follow this 
same route here, introducing a 7- or 8-dimensional parameter space to model the stellar 
bulge and bar, the stellar and gas discs and the dark matter halo. This parameter 
space is too large to be efficiently surveyed with standard scanning procedure, and in 
fact previous analysis concentrated on few sample choices of values for the different 
parameters. We overcome here this difficulty implementing instead a Bayesian approach 
to the parameter estimation (see also [16]), analogously to what is commonly done to 
estimate cosmological parameters, and also to some recent studies of supersymmetric 
models in a dark matter contest, see, e.g., [17,18]. 

The paper is structured as follows: In Section [21 we introduce the updated set of 
constraints which can be applied to the Galaxy to study the mass decomposition. The 
mass model itself is described in Section |3], while Section @] discusses the definition 
of the likelihood function we use in our analysis. Results are discussed in Section El 
concentrating in particular on the local dark matter halo density. Section [U] concludes. 

2 Dynamical constraints for the Milky Way 

Our data analysis is based on a Bayesian probabilistic inference (see [19] for a recent 
review on the subject). In a Bayesian approach to the parameter estimation, all the 
experimental informations are encoded in the likelihood function. In this section we 
introduce the data used in the construction of our likelihood function. The target of the 
analysis is the dark matter halo; datasets will be mainly constraining such component 
and are less focussed on details in the stellar or gas components. 

2.1 Terminal velocities 

We start by introducing observables which are related to the Galactic rotation curve. 
There are different methods providing information on its inner or outer portion. Circular 
velocities at galactocentric distances smaller than the Sun's can be obtained measuring 
the motion of gas clouds in the galactic plane and, under the assumption of strictly 
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circular orbits, finding along each line of sight the extreme radial velocity, i.e. the so- 
called "terminal velocity". The 2f-cm line of atomic hydrogen and the CO line have 
been used numerous times in the literature to determine terminal velocities. We take 
advantage of the modeling by Malhotra [20] of CO data from Ref. [21] and HI data from 
Refs. [22,23]. We restrain to results obtained at | sin/| > 0.35, with I being the longitude 
in Galactic coordinates, since the method assumes circular motion in an underlying 
axisymmetric potential well, and the approximation of axisymmetry is clearly not valid 
in the region of the Galactic bar, i.e. at radii smaller than, say, 3 kpc: the sample we 
will consider is then made of 50 data points at negative longitudes and of 61 at positive 
longitudes. Moreover, in fits of such inner rotation curve, we will not try to model small 
scale fluctuations, such as those most probably associated to the spiral arm structure. 
Following Ref. [13], in order to allow for random and systematic non-circular velocity, 
we assume a constant error on terminal velocities of 7 km s -1 . 

2.2 Local standard of rest velocities 

It is not possible to use the terminal velocity method to infer the outer rotation curve, 
which must be derived instead from a sample of tracers for which both distances and 
velocities are measured. We will consider the compilation of 18 regions of high- mass 
star formation, within inner, local and outer spiral arms, for which parallax and proper 
motion are measured with high precision at the VLBI, as provided in the recent work 
of Ref. [24]; we compare against the compilation of velocities with respect to the local 
standard of rest velocities, or the heliocentric velocities, versus the parallax inferred 
distances, having gauged out the estimated average value for proper motion components. 
Eleven objects out of eighteen are in the outer part of the Galaxy and will be included 
in the analysis; the rest is in the inner part and does not add information with respect 
to terminal velocities. For each object, the fit will be against the given values of the 
local standard of rest velocity ff sr , its proper motions and distance d l . 

2.3 Velocity dispersions in a tracer population 

A measurement of the rotation curve at large radii (up to about 55 kpc) has been recently 
provided in the analysis of Xue et al. [25]. As kinematic tracers, they selected ~ 2400 
Blue Horizontal- Branch (BHB) halo stars from the SDSS DR-6 with distances accurate 
to ~ 10%. Ten points in the rotation curve at radii from 7.5 kpc and 55.0 kpc have 
been derived by means of a non trivial matching of the observed radial distribution of 
radial velocities and the analogous distribution obtained from a numerical simulation of 
the Galaxy. We cannot follow this route; alternatively one can use the BHB dataset to 
trace the potential of the Galaxy at large radii $(r) (spherical symmetry is assumed in 
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all components here) using Jeans Equation for the second moment of the velocity: 



h 2-p*<T r = -P+-T- (1) 

dr r dr 

where p*(r) is the density profile of halo stars, a r their radial velocity dispersion and 
the anisotropy parameter (3 = 1 — of/of, with a t the tangential velocity dispersion. In 
case of constant non-zero /3, the solution of Eq. [T] becomes: 



1 r°° d$> 1 r°° 

—- r^ p^f)— = -—— \ dr f 2 ' 9 -V*(f)e 2 (f) (2) 

o*[r) J r dr r z P p*(r) J r 



with the circular velocity. Xue et a/, suggest to use p* oc r -3 ' 5 as determined for 
Milky Way's halo stars at 10-60 kpc from SDSS data [26]; it is also argued, based on 
their numerical simulations and for their particular survey volume, that: i) the line- 
of-sight velocity dispersion (the quantity which is measured) <J/. . s .(r) ps 0>(r); and ii) 
P = 0.37. We will assume that the first condition holds, while relax the second condition, 
assuming f3 as a free parameter. 



2.4 Local circular velocity and Galactocentric distance 

The Oort's constants, A and B, naturally appear in the expressions for the mean radial 
velocity and proper motion of astronomical objects in circular motion around the galac- 
tic center. The measure of their difference and sum constrain, respectively, the local 
normalization and slope of the rotation curve, i.e.: 

A~B = % A + B = -m , (3) 



iV ~ \dR) 



o 



were Go is the local circular velocity and Rq is the galactocentric distance of the Sun. 

Constraints on the quantity (A — B) have been derived by measuring the proper 
motions of different tracers. From the Cepheid proper motions measured by the Hip- 
parcos satellite [27] it was found A-B = (27.2 ± 0.9) kms" 1 kpc -1 . Stars of type 
O, B5 and B6 also provide good samples to determine Oort's constant because of 
their low velocity dispersion [28]. A combined analysis of the Hipparcos proper mo- 
tions and ground-based radial velocity of 0-B5 stars gave as a result [28] A — B = 
(30.06 ± 0.98) kms -1 kpc -1 . A similar analysis carried out with 0-B6 stars provided 
the value A - B = (32.0 ± 2.24) km s -1 kpc -1 [29] 0. A value of (A - B) at the high 
end side is favored also by the recent fit of the proper motions of a sample of masers 
in massive star forming regions (one of the sample quoted above for local standard of 
rest velocities) to a linear parameterization of the galactic rotation curve, which finds 
A-B = (30.3 ±0.9) km s -1 kpc -1 [24]. 



^^We combine in quadrature the errors on A and B reported in [29]. 
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The ratio Qq/Ro can also be measured independently from Oort's constants, moni- 
toring the proper motion of the central radio source Sgr A*. Assuming that Sgr A* is 
motionless and located at the Galactic center, its apparent motion is then interpreted as 
a consequence of the motion of the Sun with respect to the Galactic center. By removing 
the solar motion in longitude from the reflex of the motion of Sgr A* in longitude, the 
authors of [30] find: 

Q /Ro = A- B = 29.45 ± 0.15 km s" 1 kpc" 1 . (4) 

We take this precise estimate as our reference value for the combination A — B. 

The quantity (A + B) is less precisely determined. In the literature not even the 
sign of the slope of the rotation curve [i.e. -(A+B)] at R > Ro seems to be clearly 
understood. Positive values are quoted in [24,31,32], negative values in [25,28,33] and 
values compatible with zero in [29,34]. (A + B) has been recently estimated through 
the kinematics of old M type stars of the thin disk studied with SDSS data [35]. The 
authors find fgf = -£f=f} = -0.006 ± 0.016 which combined with Eq. ® gives 

A + B = 0.18 ±0.47 km s _1 kpc _1 . (5) 

We assume Eq. (jSJ) as our reference value for the combination A + B. 

The monitoring of stellar orbits around the central black hole has recently led to 
breaking the degeneracy between the black hole mass and value of the distance of the 
Sun from the Galactic center (setting the conversion from measured angular sizes to 
proper sizes), obtaining [36]: 

R = 8.33 ±0.35 kpc. (6) 

Combining this value with Eq. (jl]), one finds Go = (245 ± 10.4)kms _1 , i.e. the present 
observations seem to favor values of Bo larger then the IAU recommended reference 
value Go = 220 kms -1 [37]. Remind that, neglecting non circular motions, Go is the 
velocity, measured from the Galactic center, of the local system of rest, which is defined 
as the frame comoving with an observer moving along a closed orbit passing through 
Rq. Therefore, Go has to be extracted from the motion of the Sun with respect to some 
reference object, e.g. SgrA*, or a specific stellar population. We now know that the 
larger is the velocity dispersion of the reference sample of stars, the smaller is the local 
circular velocity estimated from it. This effect goes under the name of asymmetric drift. 
In principle, the correct value of Go should be therefore estimated in the limit of zero 
velocity dispersion of the reference sample of stars. For these reasons, it has been pointed 
out in [28] that the recommended value Go = 220kms~ 1 is probably underestimated 
because extracted averaging on populations of stars with different velocity dispersion. 



2 This number is extracted from a data sample analysis which assumed Oort's constant estimated 
by [27]. 
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2.5 Local surface mass density 



The total amount of matter in the solar neighborhood can be inferred by the motion of 
stars in the direction perpendicular to the Galactic plane. Kuijken & Gilmore [38] have 
examined local star velocity fields and derived a constraint for the total mean surface 
density within 1.1 kpc: 

S N <i.i kpc = (71±6)M pc- 2 . (7) 

The local surface density corresponding to the visible components only has been esti- 
mated instead with star counts [39]: 

E* = (48 ± 8)M pc~ 2 . (8) 



2.6 The total mass at large radii 

The total mass within a sphere of radius R ^> Rq can be estimated from the velocity 
distribution of the Milky Way's satellites [13]. We use measurements of the total mass 
within 50 kpc and 100 kpc to constrain our Galactic model at large radii. The values 
which we considered in our analysis are [40] 

M(< 50 kpc) = (5.4 ± 0.25) x 10 11 M , (9) 

and [13] 

M(< 100 kpc) = (7.5 ± 2.5) x 10 11 M . (10) 



3 A mass model for the Milky Way 

We present in the following our reference mass model for the Galaxy. The decomposition 
is performed describing each component as either axisymmetric or spherically symmetric. 
Depending on the assumptions, the model is fully specified by, in total, 7 or 8 parameters; 
these are listed in table CD 



3.1 The stellar disc 

The mass density of the stellar disk is sketched by a one-component thin disc defined 
by the following function: 



p d (R } z) = -^e sech 2 - with R < R dm , (11) 

£Zd \ Z dJ 

where is the disc surface density and Rd (zd) sets a scale of length in the R (z) 
direction. This form is in fair agreement with the one suggested by Freudenreich [41] 
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and fitted against COBE photometric maps of the Galaxy. Freudenreich postulated 
the presence of a hole in the inner part; such holes seem to be predicted by N-body 
simulations as a consequence of the formation of a stable bar. As none of our model 
discriminators is very sensitive to fine details in the distribution of matter towards the 
center of the galaxy, to simplify calculations of the rotation curves, we neglect such 
feature and keep in mind that, in this way, we probably overestimate the density of stars 
associated to the disc in the inner Galaxy. Moreover, Eq. ( TTTT) describes a flat disc at all 
radii and can not take into account deviations from the flat reference plane [41]. 

Four parameters are introduced in Eq. (II ip . however only two will be let vary freely. 
We fix the vertical height scale to the best fit value suggested in [41], Zd = 0.340 kpc, 
since dynamical constraints we implement are insensitive to a slight variation around 
this value; we also assume that the truncation radius scales weakly with the value of the 
local galactocentric distance R , according to the prescription R dm = 12[1 + 0.07(i? — 
8 kpc)] kpc [41]. 

3.2 The stellar bulge/bar 

The bulge/bar region is not axisymmetric and can be described by: 



p bb (x,y,z) = p bb (0) 



1.85 exp (_ s J + exp 



(12) 



where 



q 2 (x 2 + y 2 ) + z 2 



si 



(13) 



x b,Vb,Zb, lb and p bb (0) are considered as free parameters of the model. This is the form 
proposed by Zhao [42], based on COBE photometric maps, and it is in fair agreement 
with the one given by Freudenreich [41]. The second term in Eq. (TT2l) describes a triaxial 
bar while the first term represents an axisymmetric nucleus whit a power law behavior 

~*a L85 [15]. 

In the analysis, we will implement an axisymmetrized version of Eq. ffT2"j) . and as- 
sume x b ~ y b = 0.9 kpc ■ (8 kpc/i? ), z b = 0.4 kpc • (8 kpc/i? ) an d q b = 0.6. These 
(axisymmetrized) values are in agreement for the best fit obtained in Ref. [42] assuming 
Ro = 8 kpc, and then scaled to an arbitrary Rq. The choice of fixing these parameters 
is again related to the lack of observables, among those implemented, to discriminate 
among these values and small deviations around them. The normalization p bb (0), or 
equivalently the mass of the bulge/bar system at a fixed R$, is kept as real free param- 
eter. 
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3.3 The dust layer 



The distribution of the interstellar medium is assumed to be axisimmetric as well. We 
infer the average density per annulus around the galactic center of atomic and molecular 
hydrogen from the analysis of Dame [43], correcting then for helium. The vertical 
distribution is assumed to be the same as for the stellar disc. The only free parameter 
entering in the definition of this term is Ro, as we rescale again with Rq/S kpc all radial 
distances. 



3.4 The dark matter halo 

There are a few possible choices here. One possibility is to follow the scheme suggested 
from results of N-body simulations of hierarchical clustering, and sketch the Milky Way 
as the spherically symmetric radial density profile: 

Ph (r) = p'f(r/a h ) , (14) 

where, according to the numerical simulations, f(x) is the function which sets the univer- 
sal, or nearly-universal, shape of dark matter halos, while p' and ah are a mass normaliza- 
tion and a length scale, usually given in terms of the virial mass M vir and a concentration 
parameter c vir . We will adopt here the definitions: M vir = 4n/3A vir po R% ir1 with A vir 
the virial overdensity as computed in Ref . [44] , p the mean background density and R V i r 
the virial radius; and c„ r = R V i r /r_2, with r_2 the radius at which the effective logarith- 
mic slope of the profile is —2. The latest simulations favor the Einasto profile [45,46]: 



f E (x) = exp 



— (x aE - 1) 

OLE 



(15) 



with the Einasto index ranging about 0.1 to 0.25 (reference value, say, «b = 0.17). 
We will also check results for the profile originally proposed by Navarro, Frenk and 
White [47] i.e. 

W(*) = ifr ^ ? , (16) 

which is most often used. 

While the NFW profile has a 1/r singularity towards the center of the Galaxy, this 
is smoothed or partially smoothed by the Einasto index a^. Once more, our analysis is 
not focussed on the central region of the Galaxy, with the model and the implemented 
constraints that are not accurate enough at small radii. We also have a reduced discrim- 
ination power with respect to models for which the central dark matter enhancement is 
totally erased, such as the Burkert profile [48]: 

fe w = ( i + *Ki + *r (»> 



8 



This model has a density profile with a constant core rather than a dark matter cusp, 
and is phenomenologically motivated since it fits better than the NFW profile the gentle 
rise in the rotation curve at small radii which seems to be observed for many external 
galaxies, especially in case of low-mass dark-matter-dominated low surface brightness 
and dwarf galaxies [49]. Another issue which we choose not to address in the present 
version of the analysis is the fact that, in principle, one should refer to Eq. (fill) as to 
the dark matter density profile prior the baryon infall, and then model and take into 
account the feedback from the baryon settling into the center of the Galaxy. Actually 
how this happens is still a matter of debate. The simplest scenario is the one in which 
the baryon infall could have occurred adiabatically, as a smooth and slow process with 
no net transfer of angular momentum between components; in the approximation of 
a spherical system in which the local velocity distribution is unchanged as well (as 
it happens, e.g., if all particles are placed initially on circular orbits), the adiabatic 
invariants drive a rather sharp enhancement of the dark matter concentration in the 
central region of the Galaxy [50]. At the other extreme, the baryon infall could have 
happened with a sensible exchanged of angular momentum between baryons and dark 
matter, with an equilibrium configurations of the dark matter that would resemble the 
Burkert profile [51]. We choose to avoid these issues and assume that the net effect from 
baryons is at most to change the initial concentration parameter to some different final 
concentration parameter, which is the quantity relevant for dynamical observations and 
dark matter detection (this in turns weakens the link between results presented here and 
the correlation pattern between M vir and c vir seen in N-body simulations). 

4 The Likelihood function 

The Galactic model of section [3] is fully specified only when the values of its parameters 
are fixed. Determining these parameters from the data by statistical inference is essential 
to investigate halo properties such as the local dark matter density. This task requires 
a link between the data and the model parameters which, in the present analysis, is 
established by the Likelihood function. We summarize briefly below the key concepts 
and the steps that are needed to derive our results. (For a more detailed description, 
see, e.g. [52]) 

The Likelihood is the joint probability density function (pdf) to observe the whole 
set of actual data given a fully specified model. It is usually denoted by C(d\g), where d 
and g represent the data and the parameters respectively. In the statistical inference the 
parameters g can be treated either as deterministic variables or as stochastic variables. In 
the first case the aim of the data analysis is to estimate the model parameters through a 
maximum Likelihood approach. Although commonly used in the study of the Milky Way, 
this technique involves a numerical maximization of Likelihood functions and becomes 
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too involved when these are defined in a parameter space of large dimensionality. In 
the second case, instead, the analysis aims to infer probability regions in the parameter 
space of the underlying Galactic Model. It turns out that this second approach is the 
most appropriate when the Likelihood is a complicated function of a large number of 
parameters. 

Seen as a function of the model parameters, the Likelihood is not a true pdf defined 
in the parameter space. It is however mapped into a pdf for the model parameters by 
means of Bayes' theorem, which states 

where p(d) is the Bayesian evidence and in the present discussion simply plays the 
role of a normalization constant. The function p(g), instead, is the prior probability 
density function and encodes our prejudice concerning the most probable values for the 
parameters before seeing the data. The prior distribution can be crucial or irrelevant in 
the inference procedure depending on how peaked is the Likelihood and therefore on how 
many informations are carried by the used data. The pdf p(g\d) is called the posterior 
probability density function and it is the main target of every Bayesian probabilistic 
inference. It reflects the change of our prejudice about the most probable values for the 
model parameters after seeing the data. When the posterior pdf has been determined, 
then one can: 

• construct marginal posterior pdf integrating p(g\d) over portions of the parameter 
space; 

• infer posterior pdf for generic functions / of the model parameters; 

• calculate averages and variances of the desired quantities with respect to p(g\d). 

Given a n-dimensional parameter space with parameters denoted by g — gi, . . . , g n , 
the marginal posterior pdf for g 1 reads 

P(9i\d) = J dg 2 ... dg n p(g\d) . (19) 

Analogous expressions hold for portions of the parameter space of dimension larger then 
one. In the next section, for instance, we will calculate two dimensional marginal pos- 
terior pdf to look for correlations between different quantities of interest. The posterior 
pdf for any function / of the model parameters is denoted by p(f,g\d). It is related to 
the posterior pdf p(g\d) by the relation 

P(f,9\d) = 5(f(g)-f)p(g\d), (20) 
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which follows from the definition of conditional pdf. In practise, as explained in the 
appendix, / and g are discrete variables varying within a given sample. In this case, 
the delta-function in Eq. f[2"U]) is equal to one for f = f{g) and zero otherwise. Finally, 
expectation values of any function / of the model parameters with respect to p(g\d) can 
be calculated as follows 

(f(g)} = J dgf(g)p(g\d). (21) 

A subset of the functions / considered in the present paper is given in tables [2] and El 
We will call them, for obvious reasons, derived quantities. A particular emphasis will be 
attached to the case / = Pdm(-Ro) f° r its relevance in the physics of direct dark matter 
detection. Let us also mention, that in the numerical calculations, we evaluated the 
integrals introduced in this section using Markov Chain Monte Carlo methods. Further 
details on this approach can be found in the Appendix. 

We conclude this section giving an explicit expression for the Likelihood function 
used in our analysis. Different datasets will be considered as statistically independent, 
which means that the total Likelihood is given by the product of the single Likelihoods 
obtained from each set independently. We start by introducing the 7-dimensional array: 

D j = (A-B,A + B, v c (Ro), £ w <i.i kpc , £*, M{< 50kpc), M(< lOOkpc)) (22) 

which has in the entries the experimental data of section To each entry of D J one can 
associate a theoretical prediction T- 7 and an error as explained in the previous sections. 
The terminal velocities v\ er) the velocity dispersions in a tracer of population a* 1 , and 
finally, the proper motions {p,\, /j^), distances d l and local standard of rest velocities of 
high mass star forming regions v\ ST , will be treated separately (here i = 1, . . . ,m with 
m = 111 for the terminal velocities data, 9 for the velocity dispersions data and 11 
for the high mass star forming regions data). Errors associated to these datasets are 
denoted by a\ er1 a\ y a l d , a\, expand a\ sv1 respectively. If a given quantity carries a bar {e.g. 
t>t er ), it denotes the experimental value of the corresponding quantity. The same symbol 
without a bar {e.g. v[ eT ) represents the theoretical prediction for such a quantity. With 
this notation we can now write the total Likelihood as follows 

i ^ (t* - piy , i i ^ « r - 4j 2 , 1 1 ^ - K'f 

2^ + 2111^- ail 29^ 
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(23) 
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For data referring to the high mass star forming regions, we have introduced a mini- 
mization procedure since velocities and distances have high correlated errors and values 
of R are not univocally defined (see the discussion in [13]). 

Concerning the choice of the prior pdf, we assign flat priors to all the model param- 
eters. We can therefore write p(g) = Ylk=iPk(.9k) where the single Pk(gk) read 

{const for g^ in < g k < gf** 
(24) 
otherwise 

The ranges (g™ m , fi 1 ™ 3 '*) for the parameters of the underlying Galactic model are listed 
in table [TJ 



5 Results 

We can now present the results of our analysis, namely the marginal posterior pdf, 
means, standard deviations and confidence intervals for different quantities of interest. 
As already stressed in the previous section, a particular emphasis will be given to the 
local dark matter density, since it is a very important quantity for direct dark matter 
detection experiments. 

To reach this goal, we run a dedicated FORTRAN code which is partially based on 
Superbayes [53] and cosmomc [54], from which we borrowed Markov Chain Monte Carlo 
routines. 

In the analysis we consider separately different functional forms for the dark matter 
halo density; we start with the case of the Einasto and NFW profiles. The latter has 
been the most commonly used so far in the literature, while the former seems to be 
favored from the most recent N-body simulations. As we will show, however, results are 
essentially independent from the choice of the dark matter profile, since as we already 
stressed the two profiles mainly differ in the central portion of the Galaxy where our set 
of dynamical constraints are less efficient. 

Fig. (TjQ) shows the pdf of the Galactic model parameters when an Einasto profile 
is assumed. Except for the central bulge/bar energy density pb&(0) and the central 
surface mass density of the disc S^, which are poorly constrained, the other parameters 
exhibit a clear peaked distribution. In particular the virialization mass M vir is very 
well constrained within our dataset. Analogously, Fig. §2§ shows the pdf for the model 
parameters in case of a NFW profile. There is a clear peaked distribution for all the 
dark matter related parameters as well as now for both disc parameters. Fig. ([3]) shows 
marginal posterior pdf for a few derived quantities which are some of the quantities we 
implemented as inputs in the analysis. The figure clearly shows consistency between 
input and output in our procedure. Note in particular that the local circular velocity, 
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Figure 1: Marginal posterior pdf of the Galactic model parameters (Einasto profile). In each plot the 
green vertical line represents the mean of the corresponding parameter while the red dashed line its 
best fit value. The yellow (blue) bar above the curves indicates the largest interval including the 68% 
(95%) of the total probability. 
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Figure 2: Marginal posterior pdf of the Galactic model parameters (NFW profile). In each plot the 
green vertical line represents the mean of the corresponding parameter while the red dashed line its 
best fit value. The yellow (blue) bar above the curves indicates the largest interval including the 68% 
(95%) of the total probability. 
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Figure 3: Marginal posterior pdf of a few derived quantities (Einasto profile). In each plot the green 
vertical line represents the mean of the corresponding parameter while the red dashed line its best fit 
value. The yellow (blue) bar above the curves indicates the largest interval including the 68% (95%) of 
the total probability. 



which is a quantity entering critically in direct detection experiments, is rather sharply 
peaked around its central value of 243 km s -1 . 

From Figs. (T4j) and (jSJ), one can see the correlations between the parameters defining, 
respectively, the Einasto and NFW halo, as well as the local galctocentric distance. In 
the first case, the Einasto index is poorly constrained, reflecting again our lack of 
information at small Galactic radii. On the other hand viral mass and concentration 
parameter are much more efficiently constrained. The mean value of the concentration 
parameter is about 18 which stands at the upper end of the value that is predicted by N- 
body simulation for a dark matter halo of virial mass around 1O 12 M ; note however the 
possible mismatch between our result and N-body simulation results eventually related 
to the baryon infall effect. One can also see that the concentration parameter is fairly 
correlated to the value of the local galactocentric distance which in our distribution 
has a mean of about 8.2 kpc. Analogously the result for the NFW case gives a very 
accurate determination of the virial mass and a sharp correlation between virial mass 
and concentration parameter. The concentration parameter is again large with a mean 
equal to 19.7. 

Fig. (El) contains the main result in this analysis. We show the marginal posterior 
pdf for the local dark matter density pdm(Rq) for the Einasto profile (top left panel) and 
the NFW profile (bottom left panel). In the same plot we also show the result obtained 
when modelling the dark matter halo with the Burkert profile (bottom right panel), in 
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Figure 4: Two dimensional marginal posterior pdf in the planes spanned by the combinations of the 
Galactic model parameters which determine the dark matter halo in the case of an Einasto profile. The 
normalization is such that at the maximum the posterior pdf is equal to one. The black dots correspond 
to the means of the plotted posterior pdf. One and two sigma contours are also shown. 
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Figure 5: Two dimensional marginal posterior pdf in the planes spanned by the combinations of the 
Galactic model parameters which determine the dark matter halo in the case of a NFW profile. The 
normalization is such that at the maximum the posterior pdf is equal to one. The black dots correspond 
to the means of the plotted posterior pdf. One and two sigma contours are also shown. 
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Figure 6: Marginal posterior pdf for the local Dark Matter density. Top left panel: assuming an Einasto 
profile and applying all the constraints. Top right panel: assuming an Einasto profile and applying 
different subsets of constraints. Global constraints include M(< 50kpc), M(< lOOkpc) and S| z i<i .ikpc- 
Tracers constraints include the local standard of rest data, the terminal velocities and data referring to 
the high mass star forming regions. Bottom left panel: assuming a NFW profile and applying all the 
constraints. Bottom right panel: assuming a Burkert profile and applying all the constraints. Curves 
and bars have the same meaning as in the previous plots. 
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analogy to the analysis just described for the other two profiles. The result is essentially 
halo profile independent, with the 68% confidence level associated to this average only 
slightly changing in the three cases. For the Einasto and NFW profiles Pdm(Ro) = 
0.385 ± 0.027 GeV cm" 3 and p DM {Ro) = 0.389 ± 0.025 GeV cm" 3 , respectively, while in 
the Burkert case we find Pdm(Ro) = 0.409 ± 0.029 GeV cm -3 . We can therefore give a 
profile independent determination of Pdm(Ro) which has an accuracy of roughly the 7%. 

Figs. (GO) and (IE]) show two dimensional marginal posterior pdf for pdm(Ro) against 
all the parameters in the Milky Way mass model and a set of input observational quan- 
tities, in case of the Einasto profile. One can deduce that the very precise determination 
of Pdm{Rq) is not stemming from any single particular observable or from an ad hoc 
choice of the parameters defining the mass model since, apart from the obvious corre- 
lation with the local mass observables such as the local stellar disc surface density and 
the local total integrated surface density, other clear correlations are more difficult to 
single out. On the contrary, the small error on Pdm{Ro) comes from a combination of 
complementary constraints. To understand which pieces of information are needed to 
determine Pdm(Rq), it is useful to write the local dark matter density as follows 



where K includes a contribution to Pdm(Rq) from the mass density of the baryons plus 
a correction due to the fact that such a component can not be treated as spherically 
symmetric. The local amount of baryonic matter is essentially constrained by the local 
disc surface density £*. The first term in Eq. (I25l) . instead, depends from several aspects 
of the Galactic model, namely the Sun's galactocentric distance, the local values of the 
rotation curve and its local slope. Note that all this informations are mainly associated 
to local observables. We find that a combination of the precise estimates for ©o and 
A — B quoted ins section [2, together with the informations on the slope of the rotation 
curve coming from A — B, the terminal velocities and other tracers, are efficient datasets 



such as M(< 50kpc), M(< lOOkpc) and £| 2 |<i.ikpc - or constraints associated to different 
tracer populations alone, could have not lead to a determination of the local dark matter 
density as precise as the one we quote in this paper (see top right panel in Fig. (Q). In 
Fig. @ we plot, for the Einasto profile, the rotation curve calculated for the point in 
the parameter space corresponding to the mean of the Einasto posterior pdf. 

All the averages with respect to the different marginal pdf and the associated confi- 
dence intervals for the corresponding Galactic model parameters can be found in tables 
[2] and [31 In all the one-dimensional plots discussed above, the red dashed lines represent 
the best fit value of the parameter considered in that plot, while the solid lines its aver- 
age with respect to the plotted pdf. We say that a point in parameter space gives the 




(25) 




It should be now clear that global constraints alone - 
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Figure 7: Two dimensional marginal posterior pdf in the planes spanned by the local dark matter 
density and one of the Galactic model parameters in the case of an Einasto profile. The normalization 
is such that at the maximum the posterior pdf is equal to one. The black dots correspond to the means 
of the plotted posterior pdf. One and two sigma contours are also shown. 
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Figure 8: Two dimensional marginal posterior pdf in the planes spanned by the local dark matter 
density and one of the derived quantities in the case of an Einasto profile. The normalization is such 
that at the maximum the posterior pdf is equal to one. The black dots correspond to the means of the 
plotted posterior pdf. One and two sigma contours are also shown. 
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Figure 9: Top panel: Galactic rotation curve in the case of an Einasto profile. Different curves are 
associated to the contributions of the various Galactic components. 



best fit to the data when it minimizes the effective xls defined as follows 



Xeff 



-2\np(d\g). 



(26) 



Means and best fit values differ considerably when the Likelihood surface is character- 
ized by a narrow region in parameter space with a very high Likelihood surrounded by 
other regions with a slightly lower Likelihood but whose volume in parameter space is 
much larger. In the present analysis, however, means and best fit values are rather close 
in all the plots, with the exception of the already mentioned bulge/bar central energy 
density. This suggests that our Likelihood contributes to the posterior pdf more then 
the prior distribution (see Eq. (118p ). This indication is confirmed by Fig. ( jTUi) which 
shows the marginal posterior pdf for a few derived quantities obtained through a M CMC 
scan performed without imposing any observational constraint. Such pdf represent the 
true priors which are actually used for the derived quantities. Since the derived quan- 
tities depend non trivially on the model parameters, the corresponding priors are not 
necessarily flat, even if flat priors in the parameter space are assumed [18]. In our case, 
however, in the range where the Likelihoods are different from zero, it is always possible 
to disentangle the peaks of the Likelihoods from the priors on the derived quantities. 
Therefore, the observational constraints considered in our analysis are stringent enough 
to extract reliable informations from the data. 
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Figure 10: The blue lines represent marginal posterior pdf for a few derived quantities obtained by a 
MCMC scan with no observational constraints. The green dashed curve are the corresponding Likeli- 
hoods. An Einasto profile has been assumed. 



6 Conclusions 

We have produced a novel study on the problem of constructing mass models for the 
Milky Way, concentrating on features regarding the dark matter halo component. We 
have implemented a variegated sample of dynamical observables for the Galaxy, includ- 
ing several results which have appeared recently. We have also developed our analysis 
introducing a rather general scheme to describe the different mass components for the 
Milky Way, introducing a model with a large number of parameters (8 or 7 in total). 
Compared to previous studies of this kind, which did concentrate on few sample choices 
of values for the different parameters, we have studied the full parameter space by im- 
plementing a Bayesian approach to the parameter estimation, analogously to what is 
commonly done to estimate cosmological parameters, and a Markov Chain Monte Carlo 
to explore it. Our results show that this novel approach is fully successful. 

The main result of this analysis is a novel determination of the local dark matter 
halo density which, assuming spherical symmetry and either an Einasto or an NFW 
density profile, is found to be around 0.39 GeV cm -3 with a 1-cr error bar of about 7%; 
more precisely we find a Pdm(Rq) — 0.385 ± 0.027 GeV cm" 3 for the Einasto profile and 
Pdm(Ro) = 0.389 ± 0.025 GeV cm" 3 for the NFW. This is in contrast to the standard 
assumption that Pdm{Ro) is about 0.3 GeV cm -3 with an uncertainty of a factor of 2 to 3. 
Our results indicate that this accuracy is preserved even considering other spherical dark 
matter profiles such as the cored Burkert profile for which we find Pdm{Rq) — 0.409 ± 
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0.029 GeV cm -3 . Sketching how solid this result is when introducing axisymmetric or 
triaxial dark matter profiles will be the goal of a dedicated study. 

A very precise determination of the local halo density is very important for interpret- 
ing direct dark matter detection experiments (see for instance the recent analysis [55]). 
Indeed the results we produced, together with the recent accurate determination of the 
circular velocity, should be very useful to considerably narrow astrophysical uncertainties 
on direct dark matter detection. 
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A Markov Chain Monte Carlo 

To calculate the integrals which define the means and the marginal pdf introduced in 
section H] is in principle a formidable task when the parameter space is of high dimension- 
ality and the Likelihood surface exhibits a complicated structure. Monte Carlo methods 
based on a Makrkov Chain sampling are an usefull tool to overcome this complication. 
In this appendix we review this idea and show how it applies to the present analysis (see 
also e.g., [56]). 

A.l Evaluation of expectation values and marginal pdf 

Let us introduce a generic n-dimensinal parameter space which can be though of as 
the parameter space of the Galactic models investigated in this paper. In the NFW 
case n = 7, while in the case on an Einasto profile n = 8. Let us now denote a point 
in such a parameter space by g = (gi,...,g n ). Given a sample of dimension N 
(t — 0, . . . , N — 1) drawn from an opportune posterior pdf, e. g. p(g\d), a Monte Carlo 
estimate of an expectation value such as the one in Eq. ( !2~T1) is given by 



An analogous expression with f(g) = 1 and the integration voulume reduced to an 
opportune subspace leads to the desired marginal pdf. 




7V-1 



(27) 
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Thus, Monte Carlo methods reduce complicated integrals to much easier sums. The 
price to pay is the need of a smapling method. Sampling methods to obtain the sequence 
are typically based on the generation of Markov chains. Formally a Markov Chain 
is a sequence of random variables such that [56] 

p(£(-+%M, {g® :teS})= P {g {m+1) \g^) , (28) 

where S is any subset of {0, . . . , m — 1}. In other words, the first m — 1 points of the 
chain only indirectly influence the probability of g( m+1 \ In practice a Markov chain 
is assigned giving the probability p(t? in ) that a point g in is the starting point of the 
chain, and the so called transition probability T m (g,g'), which is the probability to 
find g' at the position m + 1 in the chain, starting from g at the position m. With 
these definitions, the probability of finding g* at the position m + 1 is also given by 

p m +i(g*) = T.- a Pm{g)T m (g,g*). 

One of the nice feature of Markov chains is that they admit invariant distributions, 
that is, distributions such that once they are realized at the position m* in the chain, 
they are not modified by further applications of the transition probability in all the 
subsequent positions m > m*. In general, an invariant distribution 1(g) is defined as 
follows 

I(g) = Y,I(9)T m (9,9) Vm. (29) 
a 

A sufficient condition to ensure that a given pdf Tt(g) is an invariant distributions of 
a Markov chain generated by a given transition probability is the condition of detailed 
balance, namely 

n(g)T m (g,g') = n(g')T m (g',g) Vm. (30) 

Once a Markov chain is correctly generated one can use it to estimate various expec- 
tation values by means of Eq. (|27|) . Clearly, in the present analysis we are interested in 
generating Markov chains which have the posterior pdf (|18|) as invariant distribution. 

Our FORTRAN code uses the Markov Chain Monte Carlo routines of Superbayes [53] 
which generate Markov chains by means of the Metropolis algorithm. Given a starting 
point g m drawn from an assigned distribution Pm(^in)) the Metropolis algorithm returns 
the subsequent point g out preceding as follows. First, it proposes a new point g which 
is drawn from a given distribution S(gi n ,g), called the proposal distribution. Then, the 
proposed point g is accepted with probability 

A( 9 „,S) = ,„ m (l,^). (31) 

If the proposed point g is accepted, the algorithm sets g out = g. Otherwise, if rejected, 
the algorithm imposes g ont = g\ n and proposes a new point. The number of proposed 
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points g needed to move from g m to g out is the multiplicity assigned to g m . The sum 
(1271) counts each element of the chain, say with its multiplicity. It can be proven 
that after a number of iterations large enough, the Metropolis algorithm starts sampling 
from the desired pdf. 

Concerning our choice of the proposal distribution, we followed the prescription of 
[54]. 



A. 2 Convegence of Markov chains 

A set of Markov chains can be used to infer predictions only if they have converged 
to the desired invariant distribution, which for the present discussion corresponds to 
Eq. fll8p . We outline here the convergence criterion which we adopted in our analysis, 
refering to [57] for a more detailed treatement of the subject. 

For each halo type considered in this paper, we generated 24 chains containing 50000 
points of the corresponding parameter space. Operationally, such chains have converged 
when inferences separately derived from any single chain lead to the same conclusions. 
Given I chains with q elements sampled assuming a given dark matter profile, this idea 
can be formalized as follows. Let us focus on a generic model parameter gk- We denote 
by g % £ the value associated to gk after j iterations within the i-th chain, where % = 1, . . . , I 
and j = 1, . . . , q. The within-chains variance W and between-chains variance B/q of g^ 
are defined by the relations 



B/q = J2(gi - m 

i=l 



i=l j=l 

(32) 



where g\ = i £* =1 g\ ] and g k - = i £)L=i £?=i g\°. 
The weighted average of these quantities 

^-V^lw+B (33) 

q q 

provides an estimate of the true variance o 2 gk of the parameter g^. Indeed, if all the 
chains start from points well separed in parameter space, Eq. (1331) overestimates cr^. 
On the other hand, the within-chains variance W alone would underestimate a 2 gk . This 
argument suggests to monitor the ratio 



_ a\ + B/lq 

Rk ~ — w — ' (34) 
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considererd as a function of the iteration number of the Metropolis algorithm. The term 
B/lq in the numerator is a correction to the between- chains variance which accounts for 
the sampling variability [58] . By construction, the ratio is initially larger then one. 
However, when the chains have converged and the Metropolis algorithm starts sampling 
from the desired posterior pdf, Rk stabilizes around a constant value close to one. This is 
the case of the Markov chains generated in our analysis (see Fig. [TT1) . We also considered 
a convergence criterion based on the so called multivariate scale reduction factor R. This 
approach extends the present method, which applies to each parameter g k separately, 
to the case in which all the parameters are treated together. We refer to [58] for more 
details on the multivariate approach. 
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Figure 11: Top panel: Assuming an Einasto profile, we plot Eq.(|34p versus the iteraton number in the 
Metropolis algorithm for each parameter of the model. We also plot the multivariate scale reduction 
factor [58] . Bottom panel: As in the top panel but for the NFW case. 
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Galactic components 


Parameters 


Units 


min 
9k 


„max 
9k 


Disc 




M Q pc~ 2 


400 


2000 


Disc 


Rd 


kpc 


1.2 


4 


Dust layer 


none 








Bulge/bar 


p bb (0) 


M Q pc~ 3 


0.1 


3 


Halo 




1O 12 M 


0.1 


10 


Halo 


C v i r 




5 


25 


Einasto halo 


OLE 




0.1 


0.4 


All components 


P 




-1 


0.7 


All components 


Rq 


kpc 


6.5 


9.5 



Tabic 1: Parameters of the model. See the text for their physical interpretation. 
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Parameters 


mean 


a 


low 68.00% 


up 68.00% 


low 95.00 % 


up 95.00 % 


S d [M pc- 2 ] 


1154.14 


427.43 


683.14 


1662.94 


476.17 


1943.67 


Pbh (0) [M pc- 3 ] 


1.37 


0.80 


0.48 


2.32 


0.16 


2.9 


Rd [kpc] 


2.45 


0.21 


2.24 


2.65 


2.07 


2.90 


Ro [kpc] 


8.25 


0.29 


7.97 


8.54 


7.66 


8.81 


M mr [10 12 M ] 


1.39 


0.33 


1.07 


1.74 


0.93 


2.14 




18. Ul 


o o o 
3.32 


14.51 


21.76 


12.31 


24.36 


Q-E 


V.ZZ 


U.U r 


U. 14 


n on 
U.zy 


n i i 
U.ll 


U.OO 


a 
P 


-u.zy 


U.z4 


n cq 
-U.Oo 


-U.UD 


n en 


U.lo 


Derived quantities 


mean 


(J 


low 68.00% 


up 68.00% 


low 95.00 % 


up 95.00 % 


A - B [km s- 1 kpc- 1 } 


29.44 


0.15 


29.29 


29.59 


29.15 


29.74 


A + B [km s- 1 kpc" 1 ] 


0.07 


0.45 


-0.38 


0.51 


-0.82 


0.94 


v c (R Q ) [km s- 1 ] 


243.03 


8.49 


234.68 


251.28 


225.61 


259.13 


£* [M pc" 2 ] 


46.51 


5.47 


41.05 


51.96 


35.76 


57.23 


S| z |<l.lkpc [ M PC" 2 ] 


72.16 


4.24 


67.93 


76.37 


63.87 


80.47 


M(< 50kpc) [10 11 M ] 


5.36 


0.24 


5.13 


5.60 


4.90 


5.83 


M(< lOOkpc) [10 11 M ] 


8.59 


0.64 


7.94 


9.23 


7.37 


9.86 


Pdm(Ro) [ GeV cm- 3 ] 


0.386 


0.027 


0.359 


0.413 


0.333 


0.439 



Table 2: Means, standard deviations and confidence intervals for the model parameters and the derived 
quantities in the Einasto case. 



Parameters 


mean 


a 


low 68.00% 


up 68.00% 


low 95.00 % 


up 95.00 % 


S d [M pc- 2 ] 


1042.62 


188.91 


849.17 


1231.13 


678.60 


1404.01 


Pbh (0) [M pc- 3 ] 


1.31 


0.79 


0.44 


2.25 


0.15 


2.85 


Rd [kpc] 


2.41 


0.17 


2.23 


2.58 


2.08 


2.75 


Ro [kpc] 


8.28 


0.29 


8.00 


8.55 


7.67 


8.81 


M mr [10 12 M ] 


1.49 


0.17 


1.33 


1.64 


1.23 


1.86 




19.70 


2.92 


16.59 


22.90 


13.93 


24.60 


a 


-0.30 


0.23 


-0.53 


-0.072 


-0.79 


0.11 


Derived quantities 


mean 


a 


low 68.00% 


up 68.00% 


low 95.00 % 


up 95.00 % 


A - B [km s- 1 kpc" 1 ] 


29.44 


0.15 


0.29.29 


0.29.59 


29.15 


29.74 


A + B [km s- 1 kpc" 1 ] 


0.073 


0.44 


-0.37 


0.51 


-0.79 


0.94 


v c (Ro) [km s- 1 ] 


243.75 


8.34 


235.58 


251.79 


226.11 


259.05 


E, [M Q pc- 2 ] 


46.24 


5.38 


40.87 


51.60 


35.77 


56.87 


S| z |<i.ik pc [M© pc" 2 ] 


72.13 


4.18 


67.95 


76.29 


63.93 


80.31 


M(< 50kpc) [10 11 M ] 


5.35 


0.24 


5.11 


5.59 


4.88 


5.82 


M(< lOOkpc) [10 11 M ] 


8.56 


0.53 


8.035 


9.08 


7.59 


9.65 


Pdm(-Ro) [GeV cm" 3 ] 


0.389 


0.025 


0.365 


0.414 


0.338 


0.435 



Table 3: Means, standard deviations and confidence intervals for the model parameters and the derived 
quantities in the NFW case. 
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